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Abstract. We study the nonequilibrium dynamics of a spinful single-orbital quantum 
dot with an incorporated quantum mechanical spin- 1/2 magnetic impurity. Due to the 
spin degeneracy, double occupancy is allowed and Coulomb interaction together with 
the exchange coupling of the magnetic impurity influence the dynamics. By extending 
the iterative summation of real time path integrals (ISPI) to this coupled system, we 
monitor the time-dependent nonequilibrium current and the impurity spin polarization 
to determine features of the time-dependent nonequilibrium dynamics. We especially 
focus on the deep quantum regime, where all time- and energy scales are of the same 
order of magnitude and no small parameter is available. We observe a significant 
influence of the nonequilibrium decay of the impurity spin polarisation both in presence 
and in absence of Coulomb interaction. The exponential relaxation is faster for larger 
bias voltages, electron impurity interactions and temperatures. We show that the exact 
relaxation rate deviates from the corresponding perturbative result. In addition, we 
study in detail the impurity's back action on the charge current and find a reduction of 
the stationary current for increasing coupling to the impurity. Moreover, our approach 
allows to systematically distinguish mean-field Coulomb and impurity effects from the 
influence of quantum fluctuations and flip-flop scattering, respectively. In fact, we find 
a local maximum of the current for a finite Coulomb interaction due to the presence 
of the impurity. 
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1. Introduction 

Diluted magnetic semiconductors [1, 2, 3, 4] are an important class of materials 
for the spintronics community since they combine properties of (ferro-)magnets and 
semiconductors. Thus they allow for the all-electrical control of the magnetic degrees 
of freedom of spintronic devices. The magnetisation of semiconducting devices is 
mainly caused by the interaction of magnetic impurities in the sample with itinerant 
charge carriers. To understand the microscopic details of these magnetic materials, 
reduced model systems are required in which the individual constituents are well under 
control. This allows to study fundamental questions concerning the interplay of coherent 
quantum dynamics and dissipation under nonequilibrium conditions. An ideal candidate 
for such a model system is a small quantum dot which connects metallic leads and also 
carries a magnetic degree of freedom described by a quantum mechanical spin {magnetic 
Anderson model). 

Magnetic quantum dots have been studied experimentally in ensembles which are 
particularly suited for the investigation by laser and electromagnetic fields [5, 6, 7, 8, 9, 
10, 11]. They are designed with standard lithographic methods and are technologically 
well established. Moreover, embedding individual Mn ions into quantum dots and 
studying the electrical properties is possible. [12, 13, 14]. Small quantum dots with 
few charge carriers and a single magnetic impurity may become important candidates 
for efficient high density spintronic devices. Although Mn ions in quantum dots have 
usually spin larger than S = 1/2, there are good reasons to study the low spin case first. 
From a methodological point of view, this simple model is an ideal basis to develop 
the numerical tools necessary to treat the real-time dynamics of more general coupled 
electron-impurity spin systems as well. In addition, the magnetic Anderson model serves 
as a generic model to study electronic transport through magnetic atoms or molecules 
placed between the tip of a scanning tunnelling microscope and a substrate [15, 16, 17]. 
Likewise, our model is a phenomenological basis for microscopic studies of molecular 
junctions based on organic radicals [18]. We mention the switching of the spin state of 
the central iron ion in a single molecular complex in a double layer on gold by a low 
temperature scanning tunnelling microscope [19] for which our model also is applicable. 
Finally, it also mimics features of the dynamics of electrons in quantum dots that are 
subject to hyperfine interactions with the nuclei of the host material. In that respect, 
the understanding of the electron-impurity coupling and its influence on the dephasing 
and relaxation times of qubits is essential for the experiments realized in Refs. [20, 21]. 

Nonequilibrium quantum transport in strongly correlated systems continues to 
remain a challenging problem. Especially, reliable theoretical treatments prove to 
be difficult when no small parameter is present in the system, i.e., in the deep 
quantum regime. Approximate methods are often based on advanced perturbative 
treatments such as quantum master equations [22, 23, 24], which base on Markovian 
approximations and weak tunnelling coupling. The renormalisation group (RG) 
approach [25, 26, 27, 28, 29, 30] as well as the functional RG method are able to 
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capture essential nonequilibrium features [31, 32, 33, 34]. Real time density matrix 
RG's [35, 36, 37, 38] require to represent the system plus the macroscopic reservoirs 
by a large but finite lattice. Due to the possible appearance of finite size effects, the 
maximal propagation times are limited. 

Quantum Monte-Carlo (QMC) concepts are a priori numerically exact and have 
been adopted to nonequilibrium quantum transport [39, 40, 41, 42, 43]. As opposed to 
the fermionic sign problem, which shows up in equilibrium simulations of quantum 
many-body systems, the real-time (nonequilibrium) QMC weight function is itself 
highly oscillatory and causes the dynamical sign problem. Reaching the stationary 
nonequilibrium state at asymptotic times remains notoriously difficult [44, 39, 40, 
41, 42, 43]. An analytic continuation to the "complex voltage plane" via imaginary 
"Matsubara voltages" attempts to circumvent the dynamical sign problem. However, 
the postprocessed back transformation to real times is plagued by numerical instabilities 
[45, 46]. Recently, long-time and steady state results have been obtained by a QMC 
sampling of the diagrammatic corrections to the non-crossing approximation [47] , stating 
the reduction of the sign problem. 

In this work, we adopt the method of the iterative summation of path integrals 
(ISPI), see Ref. [48] to the magnetic Anderson model. This approach evaluates a real 
time path integral expression for the Keldysh partition function in a numerically exact 
manner and is particularly suitable for nonlinear transport. Recently, this scheme has 
been carefully verified by a comparison to existing approximations in the appropriate 
parameter regimes for the single-impurity Anderson model [48]. Furthermore, the 
prediction of a sustained Franck-Condon blockade in the deep quantum regime has 
been reported based on ISPI data as well [49]. In contrast to the stochastic evaluation 
of the real time path integral in the QMC approach, the ISPI scheme calculates the real 
time path-integral deterministically. Hence the sign problem is avoided. The fact that 
metallic leads at finite temperature suppress all electronic correlations exponentially 
beyond a finite memory time window is exploited by the ISPI method. This allows for 
an iterative propagation of an augmented reduced density operator of the system to 
arbitrary long times. By construction, the technique is limited to finite temperatures 
and/or finite bias voltages and not too large system sizes. Whenever numerical results 
are converged with respect to the memory time, they are numerically exact. Recently, 
Segal et al. [50] have provided an alternative formulation of this approach in terms of 
Feynman-Vernon-like influence functionals. 

We theoretically investigate the real time dynamics of a single-level quantum 
dot containing a magnetic impurity with spin-1/2. Exchange interaction with on-dot 
electrons results in dissipation, induced by the metallic leads to the localised impurity 
in the quantum dot. The deep quantum regime, characterised by the absence of a small 
parameter is in the focus of the present article, i.e., all interaction strengths are of the 
same order of magnitude. In particular, we are interested in the nonlinear transport 
regime, where a finite bias voltage is applied and linear response theory is no longer 
applicable. The real time relaxation of the impurity spin as a function of various system 
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parameters is investigated. Due to the additional degree of freedom of the magnetic 
impurity, an extension of the ISPI scheme [48] is necessary. Although the inclusion 
of a magnetic impurity affects the single particle dynamics in the first place, coupling 
of the impurity to the electronic density on the dot renders this extension nontrivial. 
Our accurate results show that in the considered cross-over regime, the nonequilibrium 
charge current significantly influences the quantum relaxation dynamics of the impurity. 
Likewise, the back action of the impurity dynamics on the nonequilibrium current 
becomes significant. Most importantly, we clearly show that this crossover regime is 
not accessible by perturbative means. 

The structure of the article is as follows. After introducing the model in Sec. 2, 
we express the Keldysh partition function as a real time path integral in Sec. 3 and 
show how to evaluate this path sum by a deterministic iteration scheme. We calculate 
in Sec. 4 the stationary transport current and the impurity spin dynamics. Some of the 
results are compared to outcomes of perturbative methods for the appropriate parameter 
combinations. We analyse the influence of the nonequilibrium current on the impurity 
relaxation. Furthermore, we present results of the stationary tunnelling current in the 
deep quantum regime, where rate equation results are not reliable. The dependence of 
the relaxation rate on various model parameters is presented. 

2. Model system 

We extend the single-impurity Anderson model for the electronic degree of freedom 
of a quantum dot (QD) coupled to two metallic leads (L/R) via tunnelling barriers 
in order to study magnetic quantum dots, see Fig. 1 for a sketch. We assume equal 
tunnelling barriers at the left and right side, the generalisation to the asymmetric case 
is straightforward. The total Hamiltonian is given by H = H^ot + -Pleads + Ht where we 
use units such that H — l. The Hamiltonian of the magnetic dot 

#dot = Hf ot + H imp + if int (1) 

includes the electronic part Hf ot and the part if; mp to describe the spatially fixed 
magnetic impurity as well as the coupling term between electron and impurity denoted 
as H int . 

We write the electronic part of the QD as 

^dot = #d°ot + Hl t = e*4d a + Ud\d t d\d± , (2) 

<j 

where the operator d a annihilates a dot electron with spin a. The dot is formed by 
a single spin degenerate level, which can be controlled by a gate electrode that shifts 
the electrostatic potential of the dot Hence, the electronic degree of freedom can 
assume four values {0,t 5 4, d}, indicating whether the dot is empty (0), contains one 
electron with spin a =t, 1= ±1 and energy e a = $£> + a" A/2, or is in a spin singlet state 
with double occupation (d). The Coulomb repulsion is modelled via U > when the 
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(a) V 

Figure 1. (a) Two metallic leads are coupled to a quantum dot (QD) via tunnelling 
barriers. The gate electrode allows to tune the electrostatic potential $d and the bias 
voltage V induces an electron current. The QD has a single orbital electronic level 
which is empty or occupied either by one electron in spin up or down state (b) or by 
two electrons in a Coulomb-interacting singlet state (c). Moreover, the QD carries a 
localised magnetic impurity M with spin 1/2. In (b), a single dot electron interacts 
with M via an exchange interaction J, while in (c) only Coulomb interaction is present 
for double occupation. 



dot is in the doubly occupied state d. The Zeeman level splitting A might be present 
due to possible external or internal crystallographic magnetic fields. 

The magnetic impurity (quantum spin) is included via the Hamiltonian 

Hi mp = 2~^~ Tz ' ^ 

with the Zeeman energy Aj mp . Spin operators of the impurity are given in terms of 
the Pauli matrices T X)y z with r± = r x ± ir y . The impurity spin r and the dot electron 
spins o are coupled by the exchange interaction of strength J which is captured by the 
Hamiltonian 



#int = Jr z {d\d t - djdj + -(r + d|d t + r_d\d{) . (4) 



While the first term shifts the single-particle energies, the second term induces mutual 
flips of an electron- and the impurity spin, which we denote as flip-flop processes. The 
interaction vanishes for double occupation of the dot, then two electrons form a spin- 
singlet state. 

As usual, we model the noninteracting leads by i/i ea ds = Z^kap e kc] cap Ckap, where 
Ck ap annihilates an electron with spin o and wave vector k in lead p = L, R = ±1. A 
bias voltage V is symmetrically applied between the (thermal) leads and shifts their 
electrochemical potentials such that p p = peV/2. Finally, the tunnelling Hamiltonian 
Ht = Ylkap l^Ckap + 7* c fc<jp^<7 describes the energy- and spin-independent tunnelling 
of electrons with the amplitude 7. We assume a constant density of states g(e^) around 
the Fermi energy and work in the wide-band limit. Then, the tunnelling is parametrised 
by the parameter T := 2tt |7| 2 g(ey). 
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3. The ISPI approach 

In order to determine the nonequilibrium electron current and to study the relaxation 
properties of the impurity, we generalise the approach of the iterative summation of path 
integrals (ISPI) of Ref. [48]. This is a nontrivial step and requires to incorporate the 
impurity quantum dynamics consistently within the quantum many-body formalism. 



3.1. Path integral, generating function and short-time propagator 

The ISPI approach is based on the Keldysh partition function Z = Tt{Uk[H]p(— oo)} = 
1, where U K [H\ = T K exp{— % J^Hdt} is the propagator along the Keldysh contour /C. 
The time ordering operator is T K and p{— oo) is the density matrix of the system's initial 
state [51, 52]. A generating function Z[r](t)] is constructed to calculate the expectation 
value of an Operator 0(t) via 

5\nZ[r](t)} 



(O)t 



(5) 

T]=0 



S V (t) 

with Z[rj(t)] = Ti{UK[H+if](t)0} p(—oo)}. Then, Z[rj(t)] is represented by a (fermionic) 
path integral. Throughout this work, we assume an initially (ti = — oo) empty 
quantum dot prepared in the spin-up impurity state. The full Keldysh time interval 
A t := tf — ti = NS t is discretised into time steps S t . A short-time propagator Us t is 
introduced such that it is related to the exact propagator U(t + o~t,t) = U$ t + 0(S}). 
Subsequently, a complete basis of fermionic coherent states is inserted between every 
two short-time propagators, accounting for the Coulomb and the flip-flop interactions. 

To construct a particular Us t , we separate the total Hamiltonian H into the diagonal 
part H and a nondiagonal part Hi with respect to appearing dot operators d a and g?^. 
This yields H = H + H x with 

H^H^ + H^ + Ht. (6) 

In the interaction picture, the full real-time propagator from the initial to the final time 
can be written as 

°° rtf rtf 

U(t f ,ti) = J2(- l f~l •••/ dhdt 3 ...dt N -! 

x U (t f , fjv_i)ifi . . . HMh, t^H^oih, , (7) 

where the Keldysh contour is divided into N — l pieces of free propagation by Uo(tk+i, tk) 
that are connected by iV — 2 interaction vertices —iH\dtk acting during the transition 
time dtk- Here, we define t\ := ti and ijv : = tf. When tk+i — tk + S t and in the limit of 
very small 5t, the system can either propagate freely during 5t via the free propagator 
Uo(5t) = U^tk+iitk) = Uo(tk+i — tk) or undergo at most one transition induced by Hi 
within the interval < t' < St. Hence, the short time propagator takes the form 

U 5t = Uo(St) - % [ dt'U (S t - 0#i W) = -Uo^t) (1 - mSt) : + 0{S 2 t ) , (8) 
Jo 
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where the commutator [U (t), Hi] = 0(t) is used to bring H 1 to the right of U (t). The 
remaining time integral is evaluated exactly up to order 0(Sf). A comparable error is 
obtained by normal ordering, denoted by colons in Eq. (8). 



3.2. Discrete Hubbard- Stratonovich Transformation 

Next, we treat the Coulomb term H^ ot in a similar way as in Ref. [48]. Since 
[Hq, H^ ot ] = 0, the propagator for the system with Coulomb interaction factorises into 
a free part Uo(S t ) and the interacting part exp{— iH^ ot S t }. We apply the Hubbard- 
Stratonovich (HS) transformation and obtain 



exp{-iH^ ot 5 t } = ^ ex P { 



^(n t + + i(\ St (n t - n±) 



(9) 



C=±i 

with the HS parameter A,5 t defined via 

X St 5 t = sinh" 1 ^sm[U5 t /2] + isin" 1 v / sm[U5 t /2] . (10) 

The variable £ = ±1 is interpreted as a fluctuating Ising-like spin field. Note that the 
solution given in Eq. (10) is uniquely defined as long as < US t < ir. By this step, 
the interacting problem with local-in-time Coulomb repulsions is effectively mapped 
to a noninteracting problem, at the price that the appearing spin fields interact over 
a finite time range with each other. In general, the condition S t T <C 1 is needed to 
minimise the systematic error of order 5f. If St is bounded from below, however, U 
should be adjusted in agreement with the condition in Eq. (10). The exponential 
in Eq. (9) commutes with Hq and we can write the full short-time propagator as 
exp{— i(Ho + H^ ot )S t } = 1/2 X<=±i ex P{ — * Ho$t}, thereby absorbing the classical part 
of the Coulomb interaction into the single-particle dot energies according to 

ei(S t ) = 6 a + ^ + i<j(\s t . (11) 

In passing, we note that due to the imaginary energy component, Hq should not be 
considered as a Hamiltonian. Instead, we obtain a short-time propagator by enforcing 
normal ordering, i.e., 

U l--=\H ■■eM-iH&t}: ■ (12) 
C=±i 

Combining the short-time propagators again into the full path integral, the path sum 
extends over all tuples {(} := ((jv, ■ ■ ■ °f the HS fields. 



3.3. The remaining interaction terms 

The tunnelling term H? is quadratic in the number of fermions but contains both dot and 
lead operators. Therefore, the stationary state of the isolated system is in general not 



Nonequilibrium quantum dynamics of the magnetic Anderson model 



8 



an eigenstate of the system with tunnelling. However, for arbitrary electronic coherent 
and impurity states \^/ T ) = |^)|r), the matrix elements can be written as 

(tf r '| : exp{-^} (1 - iH T 6t) : \^ T ) = (V T '\ : exp{-i(# C + H T )5 t }: \^ T ) + 0(5f) . 

(13) 

Hence, the short time propagator is obtained by adding the coupling term Ht to Hq in 
Eq. (12). To be specific, we introduce the total electronic coherent state as 

i*> = n^ 1 - wo n^ 1 - ^4)io> , (14) 

kcrp cr 

where Q a and z hav are Grassmann numbers for dot and lead electrons. 

Flip-flops of the electron- and impurity spin are mediated by the non-diagonal 
exchange coupling term H^ t in Eq. (6) . According to Eq. (8) it enters the short time 
propagator as 

U * = \Y,'- ex P{-*(#o + H T )6 t } (1 - iH&St) : . (15) 
C=±i 

In fact, this structure of the short-time propagator motivates our choice of a "mixed" 
basis of electron coherent states and impurity states |r). Most importantly, a 
straightforward derivation of a path integral in this basis becomes possible as the 
paths are separated into parts that contribute to the matrix element (\I/' r '|?7,5 t |\l/' r ) 
either for aligned (oc 1) or opposite (oc H^ t ) impurity spin orientations. This form is 
particularly useful with respect to a numerical summation over discrete impurity paths. 
An equivalent short-time propagator in form of a single exponential : exp{— iH^St}: 
with := Hq + Ht + H^ t is much more cumbersome. 

In contrast, the short-time propagator 1 — iH^5 t , though exact up to 0(5%), is 
not convenient to construct a path integral. Consider, the case of opposite impurity 
spin states r ^ r' . The phase accumulated by the free propagation between two 
instantaneous flip-flop events is missing, i.e., 

<^'|1 - iH<8 t \V)^ T , = r+i^t + <W-i*W)e*'* . (16) 

Such a propagator would only be correct if the transition process lasted the whole time 
span 5 t instead of being instantaneous. In the resulting path integral, the system does 
not propagate freely between consecutive flip-flop processes with the resulting continuous 
limit 5 t — > being unphysical. 

3.4- Constructing the full path integral 

The path integral for the generating function Z\rf[ is obtained by using Eq. (15) and 
the corresponding Grassmann fields \I/ as well as the discrete paths {r} and {(}. This 
yields 

m = E fam {-^Y p\{r}v s ^^ . (17) 
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where the path sums over impurity and HS spin-fields are performed over the 2iV-tuples 
{tj} = (j2N, ■ ■ ■ , Ti) and {Q} = (C2N, ■ ■ ■ , Ci) with Tj, Q = ±1. Within an impurity path 
{r}, m flip-flop transitions occur on the Keldysh contour, where I of them lie on the 
lower branch. We discuss the building blocks of Eq. (17) in the following. 

In Eq. (17), the action S = Si mp + S e \ + So[v] is given by the sum of the free 
impurity action and the electronic action. Due to possible source terms, depending 
on the observable of interest 0(t), S [v] is added when necessary (see below). The 
electronic action contains all coupling terms of the electrons to the impurity, to the 
leads (via the tunnelling) and to the HS fields. In particular, we have the action of the 
free impurity 

A A N A r 

S imp = -^£(r fe -r 2Ar _ fc+ = -^ J dtr(t), (18) 



k=2 



with the discrete and continuous representation given as the first and the second 
expression, respectively. This equation also illustrates how a well-defined continuous 
limit of a discrete path can be obtained (cf. Fig. 2). The electronic action S e \ can be 
represented as 



'el — Sfot + Pleads + Sj 



[ dtdt>V(t)(G%lW) (19) 

JK 

with the inverse electronic Keldysh Green's function (G el )^ naturally given in terms of 
the action S e \ with the three contributions 

Sfot = I dt* a (t)[id t -E a (t)}* a (t), 

a ^ K 

"Sleads ^ / dtc kap (t)[id t - e k }c k ap(t) , 

S T = ^ I dt[^ a (t)c kap (t) + Y~Ckap(t)Mt)} ■ (20) 

kap ^ K 

We note that one of the time integrations in Eq. (19) is trivial, since (G 01 )^, is 
proportional to S(t — t'). Since the bias voltage enters through the respective lead 
equilibrium density matrix, see Eq. (23), electronic energies are the same in both leads. 
We have used the effective dot energy in Eq. (20) defined as 

E a (t) = E+{t) = e CT + I + Jar(t) + ia((t)X St , (21) 

on the forward branch, while E~(t) = E+(t)* on the backward branch. 

The polynomial -P[{t}] in Eq. (17) depends on the impurity path {r} = {t + }({t~}) 
for the forward (backward) branch of the contour. Then, we collect all indices of the flips 
into the tuple T^ p = (k^_ e , . . . , kf) (sorted in ascending order) along the forward path 
{ r +} : = (tjv, ■ ■ ■ , Ti) with r k + ^ r fc+ _ 1 for all k + E T fl + p . Accordingly, T fl 7 p = {kj,..., k^) 
is the tuple of ascending flip indices along the backward path {t~} := (t 2 n, ■ ■ ■ , Tjv+i) 
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Figure 2. (a) Exemplary impurity path (blue line) for which the flip-flop polynomial 
is constructed. The Keldysh contour is divided into N — 1 = 8 segments of length St 
between 2N = 18 time vertices. The impurity path (tuple of black and red arrows) 
realizes m = 8 flip-flops with five flip-flops on the forward and i = 3 flip-flops on 
the backward branch. First, the flip index tuples T flip are constructed by assigning to 
each flip-flop the index of the Keldysh time that is later with respect to the real time. 
Hence, if two consecutive spins have opposite orientations, the corresponding flip-flop 
has the time index of the spin on the right-hand side of the flip (marked in red). We 
have = (8, 6, 5, 4, 2) and T^ p = (14, 12, 11). The polynomial written in (b) follows 
upon using Eq. (22). We note that the electron's spin flips in the opposite direction as 
compared to the impurity (not shown). 



with Tfc- 7^ Tfc-+i for all k" G T$ . Note that a flip index on the backward path is 
labelled according to the smaller step index of the flipping spins corresponding to the 
later time. The impurity polynomial can be expressed in terms of the Grassmann fields 

as 

p[{r}]:= If IT < 22 > 

Figure 2 illustrates an example of an impurity path. The fields of the forward and 
backward parts of the continuous inverse Green's function (G el ) _1 are not coupled since 
the corresponding Hamiltonian is diagonal. However, the associated discrete version 
connects fields from both branches via the upper right element (1,2 AT), which is due 
to the system's initial state p(ti), see Ref. [52] for details. Throughout this work, we 
assume factorising initial conditions 

p(U) = |0,Tj)(0,Ti| PlPr, (23) 

where p p oc exp{— f3 J2ka( e k — A*p) c fe CT p c fco-p} * s ^he equilibrium density matrix of lead 
p = L/R at temperature T with = (k^T) -1 and |0,Tj) denotes the empty dot with 
the impurity in the initially prepared orientation r» = | T)- 
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When constructing the action So[i]] f° r an observable 0(t m ), evaluated at the 
measurement time t m , we replace every instance of an impurity- or electron operator in 
the observable by the corresponding spin- and Grassmann fields from the forward branch 
at the time step closest to t m . This choice is arbitrary and a replacement on the backward 
branch would not change physical results. Let us assume that t m — t{ — {km — l)^t- 
The electron operators are replaced by their Grassmann field with time index k m . 
Correspondingly, we substitute the Pauli matrices according to r v (t m ) ^ 4 km) with 
ri fcm) := (1 - r fem r fem _ 1 )/2,rf m) := -i(r km - r km ^)/2 and rf m) := r km . Since the 
matrix elements (T'\T Xty \r) are non-zero only for r ^ t', the Pauli matrices should be 
replaced by field expressions that include neighbouring spins. In other words, only if 
a flip-flop occurs at time t m , the fields are non-zero. On the forward branch, a 

flip-flop with Tfc = — Tfc_i is associated with time step k. Then, So[v] '■= vO and 

(O)(t m ) = -td rj lnZ[ V ]\ rj=0 . (24) 

A generalisation to observables with two and more time parameters, e.g., correlation 
functions, is possible via higher-order derivatives of the generating function. We 
calculate the charge current I(t rn ) via the source term 

kcrp 

and the expectation value of the impurity (r z (t m )) from S Tz = i]T km . 
3. 5. Tracing Out Electron Degrees of Freedom 

Next, we perform traces over the lead degrees of freedom and perform the path integral 
over all Grassmann fields Cfc CTp , tk cp . Contour time integrations are transferred to their 
respective real time counterparts. This results in two integrals over real time for the 
(+) and (— ) branch and generates the 2 x 2-matrix structure for the Keldysh Green's 
function. The resulting generating function is written as a path integral with the effective 
electronic dot action S e \ = S^ ot + S env with the (effective) environmental action 

dt / dt'v a (t)~f( P ,t- o{i + ^ntUO - Sm(t)]Wt') ■ (26) 

Here, we have introduced the time non-local Keldysh matrix 

Y e -i/h>(t-f) ( —\ 1 \ 

^ f - i,) = ^ sinh[,( t -^] (l -I]" P7) 

This equation only holds for t — t' ^ 0. The singularity for t — t' will be addressed below. 
The exponential decay of 7 for \t — 1'\ — > 00 is the cornerstone of the ISPI method, which 
allows us to truncate certain long-time correlations (see below). Equation (26) can also 
be given in terms of an environmental Green's function 

/OO fOO 
dt dt'V a (t)(G env )lX(t') with (G env )-: = (G° nv )- t : +v( G Dlt> ( 28 ) 
-00 J —00 
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and 



. j .-i er sin[e\/(t-tQ/2] / * m (f) - * m (f ) . 
l u envJ M , 2/3sinh[7r(t-t')//3] V MO ! 



This expression is well defined for all values of t and t'. An explicit expression for 
(Genv)t,*' is given below in Eq. (31). 

Before we proceed, we address the (t — t') 1 singularity of (Gg nv ) _1 . In contrast 
to its inverse, the Green's function G® nv itself is finite at t — t'. Still, matrix elements 
decay exponentially with growing time differences. The calculation of observables does 
not suffer from this singularity either, since it cancels in the fraction 

^-^(OXU, (30) 



t)Z[0] 

which is used to numerically evaluate (O). Since the divergence does neither depend 
on the paths of the impurity- and the HS-spins nor on rj, we may rescale Z[rj\ by 
the singular factor without affecting (0)(t m ). For clarity, we keep the same notation 
for the rescaled generating function. The singularity originates from (Gg nv ) -1 . We 
collect all non-interacting contributions, which do not depend on rj into the sum 

E.(Gt)5 = EMolX} + (GLJS with 



(«=^-Op - ^v), (3D 
with := (e a + U/2). The Fourier transform of (Gq.o-)^ * s obtained as 

(<?o„W = 2 ^(- - w ) [ lV[2 _ _ w + u u + ir[i?(w) _ !] J ■ (32) 

where we have introduced F(u) = /l(^) + /r(^) as the sum of the two lead Fermi 
distributions. This 2 x 2-matrix is inverted algebraically and transformed back into 
time space by complex contour integration [53] or numerical integration. The function 
(Gf a ) ++ (t — t') is shown in Fig. 3. Whenever this divergence arises, we will remove it, 
see for instance in Eq. (42), by multiplication of — idetGf a , or equivalently, we replace 
i(Gf)' 1 with 

DM := GlAGfr 1 = 1 + + *K) , (33) 

with the self-energy 

JT k -i( k \*(5 t )) St - (34) 

The particular form of depends on the observable O. When O = I, we identify it 
with (Genv) -1 of E q- (29) and otherwise with (G )" 1 . 
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(t-f)/0 



Figure 3. (Gg 1 cr ) ++ vs. time difference i — t' for two combinations of u>^ and for 
f3T = 1/5. Except for t = t', the function is smooth for all real times. The imaginary 
part is discontinuous at t — t' , which turns out to be harmless for our numerical scheme. 
A possible way to cure it, is to introduce a high frequency cutoff cxp{— \lo\ /uj c } in Eq. 
(32), see dotted lines in the inset for the same vertical scale for w c = 100E 



dt' (G$) t , t , ^ (G e ') tkitl with tk = t x + {k-l)8 t . (35) 



The discrete form of the matrix D a follows via the relation 

^ rt k +St/2 pti+St/l 

(Go,o)ki = -~2 dt 

°t Jt k -8 t /2 Jti-St/2 

To remove the singularity at k = I, we choose the regularisation 
(Gt a )ii = [{G$ )(t- t f)-+ - + {Gt a )(t-t')-*o+)/2- 



(36) 



Alternatively, introducing a high frequency cut off exp{— \ou\ /u c } in the Green's function 
[48] yields a consistent result, see the inset of Fig. 3(b). For = 3T and eV = 4T, Eq. 
(36) yields (Gf^ij ~ 0.3384i and with the cutoff method with u c = 100r, we obtain 
0.32452 (the difference of ~ 4% decreases for larger u c ). However, using Eq. (36) has two 
advantages: first, it does not modify off-diagonal Green's matrix elements, and, second, 
we do not need an additional parameter uj c . 

We emphasize that both methods of regularisation obey the necessary causality 
relation, (Gf >a ) +/ + {G^Xf ~ ( G o!JzV ~ ( G o!J M' + = °- This follows from the causality 
structure of the Green's matrix and the self-energies £ [52]. Here, the diagonal elements 
of both Gf a and the S matrices have to be understood as the average, i.e., the integral 
of the time non-local matrix elements in an interval 5 t around the point t — t' — 0. The 
discrete version of EJJ follows accordingly. For the current, e.g., using Eq. (29) yields 



( E 9 



eT5 t sm[eV(k-l)5 t /2) 
2/3 smh[ir(k-l)6t/l3] 



3k,m Ol,m 3k,m 



(37) 



We note that also source terms for observables O containing dot fields may be 
added to the action. The effective full inverse Green's function (G )~ l is given either by 
(G^ ot )- 1 + (G , env )- 1 for the current or {Gi ab )~ 1 + i n{G )- 1 + {GlJ)- 1 for other observables. 
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Figure 4. Exemplary impurity path (blue) with N — 9 and m = 4 flip-flops along the 
Keldysh contour. 



Plugging in those pieces, the remaining path integral is recast as a discrete sum, i.e., 



S, 



off 



t /-^ J _ J —oo J —oo 

= E^twDn^^^^c}^])- 1 }. (38) 

i>,C} 

Hence, we have obtained the Keldysh partition function as a sum over expectation values 
of the polynomial P of Grassmann numbers in a system with Green's function G e ^. The 
next step is to derive an expression for (-P[{t}]) that refers to G e ^ only by applying 
Wick's theorem [52]. With Eq. (22), we have 



[P[{r}\) = ( If *?«-n n °-^">- 



(39) 



For an odd number m of flip-flops the expectation value vanishes, since each process 
contributes a creator and an annihilator for electrons with opposite spins. Odd m implies 
an odd number of alternating products of creators and annihilators, the number of is 
different from the number of 0. When applied to any state \ip) in the trace, this changes 
the particle count by 1 and the projection with vanishes. Therefore we have to 
consider paths with even m only. As an example, we evaluate Eq. (39) for the impurity 
path shown in Fig. 4 before we turn to the general formalism and arbitrary paths {r}. 
The exemplary path features m = 4 flip-flops and we find 



(P) = (5f 5J^ B ^J5>?> = <0 t 15 D?^> <^>K> 



P^ 



pi 



(40) 



S e R and the expectation value of the mixed operator product factorise with respect to 
the spin degree of freedom. Even in the presence of spin-mixing terms this factorisation 
remains valid. Applying Wick's theorem to Eq. (40) yields (P a ) = — det S CT with 



Nonequilibrium quantum dynamics of the magnetic Anderson model 



15 









^h ++ (t-f)\ 




\^ a exp{-vr(t - t') / (3} 
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Figure 5. Dependence of j ++ on t — i' for a = T/(2/3) and = An. We depict 
the absolute value and the real as well as the imaginary part. For n\At\/j3 <J 2, the 
absolute value decays exponentially with a decay time proportional to T^ 1 . The period 
of the oscillations in the real and imaginary parts is determined by V. 



(pt) = _ detl 



(Gf S )2,7 (C* ff )2,15\ j /D i\ j J (^T ff )6,3 (Gf )y | 

J 8 > 7 l (j t' ,8 > 15 / yitrj, Jl4,3 ;i4,9 y 



The general procedure to construct S CT for an arbitrary path {r} having an even number 
m of flip-flops is as follows, 

(i) Construct the flip-flop polynomial P°"[{t}]. 

(ii) Assign indices q x < . . . < q m / 2 of annihilator fields d^ 1 , . . . , d g a m/2 that appear in 
-P ct [{t}] to the rows of 

(iii) Assign indices r\ < . . . < r m / 2 of creator fields T) r J , . . . , 0^ m/2 that appear in -P ct [{t}] 
to the columns of H CT . 

Using matrix element (E a )k,i = — *(&* fc ^C) = {G e a)q k r v ^ ne & n8 ^ expression for the 
generating function follows as 

JS f 



Z[T)} = hm ("I)' (^) m exp{zS imp } ]J det i(Gf)- 1 det 5 CT , (42) 



where the summation over impurity paths is restricted to tuples {r} with T\ = t 2 n = 
i.e., correct boundary conditions along the Keldysh contour are imposed. The limit 
St — > appears explicitly here, since there is no continuous measure used for the discrete 
spin paths, neither for the HS- nor for the impurity spins. 

3. 6. Iterative summation of the path integral 

The exact generating function in Eq. (42) is intractable due to the exponentially growing 
size of the matrices for long propagation times and an adept numerical treatment is 
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necessary to proceed. The off-diagonal elements of (G efr ) _1 , given in terms of ~f{p,t — t r ) 
in Eq. (27) decay exponentially with increasing distance t — t' from the diagonal, at finite 
T and/or V [48]. This fact is illustrated in Fig. 5. In addition, the bias voltage induces 
oscillations. The correlation or memory time r c =: (K — l)S t determines the range of the 
lead-induced correlations, they are exponentially suppressed for time differences larger 
than t c . This allows us to restrict the path summation to a finite memory time window 
and thus, the number of paths that contribute to Z[rj\, originating from the magnetic 
impurity as well as from the Coulomb interaction, is drastically reduced. We use K as a 
memory time parameter in the ISPI scheme in an extrapolation procedure for r c — > oo. 
The latter eventually gives converged results independent of r c . 

To proceed, we rotate the basis unitarily, thereby rearranging the matrix elements 
of (G^)" 1 such that increasing distances with respect to the diagonal correspond to 
increasing time differences. For a given K, we obtain the block band-matrix 



(Gf) 



eff\-l 



A 


a 


A 




1,1 






<7 




A 


2,1 


A 



1,2 



2,2 



\ 



V 



A 



N C ,N C 



N C -1,N C 
a 

N C ,N C 



(43) 



where Nc = N/K and the blocks | A \ h , are K— dimensional matrices, whose entries 

are given by those of (G^. ff ) _1 taken from the rows and columns in the range of 

{(j — 1)K + 1, . . . , jK} with j = k, I. Since we neglect exponentially small components, 

the [a] -blocks in the upper (lower) secondary diagonal have an upper (lower) triangular 

structure. The E a matrices in Eq. (42) naturally inherit the same block structure with 

the modification that the corresponding blocks [~B~\ h l are in general not quadratic. Their 

dimensions are determined by the number of flip-flops within the respective time interval. 

To illustrate the scheme, a particular impurity path is shown in Fig. 6, which is 

of length 85t, consisting of 2N = 18 vertices having 12 flip-flops. The calculation of 

the determinant in our scheme needs quadratic block matrices on the diagonal, which 

we obtain again by a unitary transformation. The off-diagonal blocks are reshaped 

accordingly. In Fig. 6 this procedure is illustrated for with r c = 2S t (K = 3), obtained 

after the rearrangement. The hatched matrix elements are disregarded. The path with 

N = 9 is divided into N c = 3 segments with K vertices on each branch. In analogy 
i , r--U 



to the blocks | A\.., the diagonal blocks | B 



contain all matrix elements G? nr , with 



(i — \)K < q,r < iK (dashed boxes). Since the number of flip-flops is odd, all the 

Instead, this particular case yields 2 x 1-matrices 



B 



matrices 
Sit , and [B 



i.i 



are not quadratic. 

3,3' 



and a 2 x 4-matrix B 



i 

2,2" 



The quadratic blocks required for the 



iteration scheme are obtained by reassigning the earliest and latest creator fields 
and Dj_ i6 from the second segment to the first and third, respectively (blue arrows). Such 
a reordering is always possible and renders all diagonal blocks of S CT quadratic. 
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Figure 6. An example of an impurity path with 12 flip-flops (top) and the associated 
S^-matrix (bottom) for r c = 2St (K = 3). The discrete path (arrows on the vertices) 
has a length of 85t (N — 9) and is divided into Nc = 3 segments of length K (separated 



by dotted lines). Depending on the flip distribution, the diagonal blocks B 



general not quadratic (boxes with dashed frames) and their determinants do not exist. 
To render them quadratic, we reassign the fields of the flips closest to the segment 
borders (blue arrows). The hatched elements are exponentially small and neglected. 



We define the recursive notation X = A, B to compactify the computation of the 
determinant for the blocked Keldysh partition function as 



X 



X 



D 



with 



X 



D 









(1 




I h+l,i 1 






Xd-i 





and 



Xi 



X 



D,D ■ 



(44) 



The double line denotes matrices which themselves consists of blocks with the subscript 
D giving their dimension in blocks. The determinant of X is calculated iteratively [48] 
in D — 1 steps of which each performs the following manipulations: 

(i) Perform a Gaussian elimination of the block in the second row, first column. 

(ii) Neglect in the resulting element of second row, second column products like 



Xjk-i fcL^Jfe ifc+i' wmcn conne ct segments beyond the nearest neighbour, 
(iii) Expand the determinant after the first column, thus reducing the problem by one 
in block dimension. 



While step (i) and (iii) are exact algebraic operations, step (ii) is the second building 
block of the ISPI method, necessary for the scheme to remain consistent with neglecting 
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correlations beyond r c . A step k — > k + 1 is performed solely based on the determinant 
after step k and the spin orientations in segments k and k + 1. Here we stress that within 
the time window r c , the ISPI scheme takes into account important non-Markovian 
effects in a natural way. In the present work, these correlations are lead-induced. 
Within a typical Markovian approximation, the real-time dependence of the GF in 
Fig. 6 is replaced by 5(t — t'). Including terms in the iteration that connect segments 
beyond nearest neighbouring i^-blocks would require information about impurity spins 
in "earlier" segments < k, which are beyond r c . After step (iii), we arrive at 



det X = det X 



X 



i,i 



X 



2,2 
3,2 



X 



2,3 



X 



(45) 



D-2 



where X 



2.2 



X 



2,2 



1,1 



X | 2 JX | n 1 [Aj 1 9 . Subsequent iteration gives the final relation 



1.2 



D 



detX = detQO^ndettfZL " S.-iS- m-S^J 



(46) 



Next, we construct the E^-matrices for finite correlation times starting from Eq. 
(42) by filling the entries with the elements of G^[{r : £}]. The latter depends on the 
entire spin path {r, (} between ti and tf. In principle, inversion of the full inverse 
Green's function is possible but out of reach for practical applications since the numerical 
effort grows exponentially with propagation time. To remain consistent with the finite 
correlation time approach, we have to find approximations of | B \ i ^ ±1 , that depend on 



spins of the nearest neighbour segments i(±l). We observe that blocks on the secondary 
diagonals contribute much less than those on the main diagonal. Hence, we expand 
(G^f) -1 in powers of the off-diagonal blocks, yielding 



G 



k,k 



A 



-i 

k,k 



and G 



k,i 



k.i 



A 



k.i 



(47) 



for all 1 < k : I < Nc and \k — l\ = 1 (the index a was omitted). The blocks [G 



kd 



are 



defined in analogy to the | A \ -blocks and the | B \ -blocks are filled as described above. 

Next, we use from Eq. (33) the relation G^ = D~ 1 Gq (T , see Sec. 3.5. The S-matrices 
are free of any singularity as well and from Eq. (47), we obtain the approximate inverse 
of D a . We multiply the result with the free Green's matrix in block form. For step 
k — 1 — >■ k, we obtain 



G 



G 



G 



G 



fc-l,k-l 

k-l,k 

k,k-l 



fe-l,fe-l 



-1 



fc-i,fc-i 



k-l,k 



k,k\ G ° 



fc-i,fe} 



k,k-l 



k,k \ 



k,k-l 



fc-l,fc-l 



G 



(48) 
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Figure 7. Exemplary paths for the impurity (a) and the HS-field 

(b). For t c — 2S t (K — 3) both are decomposed into segments of length 
if = 3 in real time, which contain 2K — 6 spins. For example, path segment 
Mi = ( T 3^ r 3 + > T 2~; r 2 + : r r> r i + ) = (t,4-,4-,t,t,t)- Accordingly, we find that 
{C}3 = (C^.C^CfTfCs^Cf .C*) = (t 5 t 5 t,4',t,4')- The compact notation combines 
the impurity- and HS spins as {r,Qj, e.g., {r,C} 2 = (t, t,4-,4-, t, t, 4-, 4-, t, t, 4-,4-)- 



This allows us to construct the B_ -blocks without calculating the inverse Green's 
function explicitly. 

Collecting all parts, we can finally express Z[rj] iteratively as 



Z[r]] = 2j Z Nc where Zj = kjj_ x Zj-\ . 

{t,(}n c {T>C}i-l 



(49) 



The real time propagator of the ISPI scheme is introduced as 

A ■ 1 U n ^{SL-HL .d 



cr,-l 



X 



a X=B,D 

with the chosen initial configuration 



Z 1 = F 1 J[ J] det(T 



1,1 



(50) 



(51) 



cr X=B,D 



Furthermore, we use the definition {r, (}j = {rJ K) . . . , rjj^ 



' • '•' 1 1' Cjjf) • • • 5 C(J-l)J(-+l) 



as 



the tuple of those impurity- and HS-spins that lie in the j-th path segment of length K, 
see Fig. 7. The propagator Ajj-x (and matrix blocks) depends on all HS- and impurity 
spins in segments j — 1 and j. The prefactor 



F j = 2 



-2K, 



e im p 



(52) 



is related to the number and the position of the flip-flops, mj is the number of flip-flops 
in segment j, out of which £j lie on the backward branch. The phase is defined as 



imp 



jK 

£ 



T, 



(53) 



/=(j_l)iC + l 
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Figure 8. (a) and (b): Trotter extrapolation (solid lines) for St — > for t c = 1/r, 
t m T = 4, f)T — 1, and = A = Aj mp = 0. (a) Charge current in units of Iq = eT/h 
for J = 0, (b) mean impurity spin polarisation for eV = 0.5T and U = 0.5T. The 
impurity spin was initially in the spin- up state Tj = 1. Tiny deviations point to 
negligible unsystematic errors. For all cases in (a), the standard deviations are below 
1%, while in (b), they rise from around 1% for J = 0.3r to about 10% for J = 0.9F 
(c) and (d): Memory extrapolation r^ 1 — > for (c) the current and (d) the impurity 
orientation. For all combinations of eV and U, the standard deviations to the linear 
fit are below 1%. Within the error margin, the numerical value of the current for 
eV = 0.6T and U = coincides with the Landauer-Biittiker value Jlb ^ 0.594/ - 



3. 7. Extrapolation procedure 

By construction, the numerical value of Z[rj] contains two systematic errors, namely 
the finite discretisation time step 8t and the finite correlation time r c = (K — l)5 t . In 
the limit 8 t — > and K — > oo, however, the iterative procedure yields an exact result. 
The major benefit of this iterative scheme is that the numerical costs scale linearly with 
evolution time t m — ij. The systematic errors are eliminated [48] by an extrapolation 
of the numerical results to vanishing Trotter increment 5% and infinite memory time 

T c — > OO. 

Due to the Trotter time discretisation, all expressions are by construction exact 
up to order 5f terms. For a fixed r c and small enough values of 5t, we extrapolate 
the numerical values of some observable to St — > and thus completely eliminate the 
Trotter error. An example is shown in Fig. 8 (a) and (b) for the current and the 
impurity orientation, respectively. Depending on the observable, Trotter convergence 
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1.0 1.5 2.0 20 40 60 80 



l/(r c r) t m T 

Figure 9. (a) Mean impurity polarisation vs. t" 1 for different measurement times 
t m for $ D = -r/2, A = T/2, J = -r/2, /3r = 5 and A imp = U = 0. Notice that 
the vertical scale is different for each t m to stress the relative variations. For all plots, 
the lower bound of the scale is set to r" lm = 0.025, while the upper bound r™ ax is as 
indicated (solid lines are guides to the eye). The dashed line illustrates the results for 
tmX = 80 (green triangles) when scaled as the data for t m T = 4 (blue crosses). The 
optimal values of (r 2 ) evaluated at the local minimum are shown in (b) as function of 
t m . A fit to an exponential with a relative standard deviation of 5 • 10 -3 is represented 
by the solid line. 

may be achieved on different scales [57]. Note that one source of errors is the numerical 
derivative in Eq. (30) which results in tiny imaginary parts of observables, ranging 
between 10~ 5 to 10~ 3 . For typical parameters, it is at least one order of magnitude 
smaller than the numerical error from the linear extrapolation. 

For the memory extrapolation r c — > oo, we do not have a strict mathematical 
argument at hand, in contrast to the Trotter extrapolation. Whenever results are 
convergent, however, we find empirically two typical behaviours: (i) either the numerical 
results for (0)(r c ) depend linearly on l/r c with small deviations, (ii) or their dependence 
on 1/t c is reasonably smooth and exhibits a local extremum. The latter case is consistent 
with the principle of least dependence, see Refs. [48, 49, 54, 55] for a verification when 
compared to analytical results. An example of the linear scaling of the numerical results 
with t~ x is illustrated in Fig. 8 (c) and (d). When (0)(r c _1 ) shows a weak dependence 
on t~ 1 in a certain corner of parameter space, we still try to apply criterion (ii). Such a 
behaviour results from a trade-off between accuracy and computational costs. A minimal 
Trotter error requires minimal 5t, while, at the same time, a maximally large correlation 
time r c is desirable. Naturally, these requirements are limited by the exponentially 
increasing numerical costs. This is illustrated in Fig. 9. 
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Figure 10. Current (crosses) and impurity spin polarisation (circles) for t m T = 4, 



I3T = 1, $ 
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0, S t T = 1/3, K = 4, and increasing values of the flip-flop 
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Convergence is obtained already for m 
flip-flop number is 2(K — 1) = 6). 
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3.8. Restricting the number of flip-flops in the memory window 

In order to reduce the exponentially growing number of contributing paths (~ 4 K ) 
without affecting the accuracy, we may exploit that Fj ~ (J5t/2) mj in Eq. (52). The 
number rrij of flip-flops in path segment j is < rrij < 2(K — 1). We observe that 
the smaller the weight of each segment is, the more flip-flops it contains. On the other 
hand, the number of path segments {r}j with rrij flip-flops (given by 4Ct ( f " 1} with 

= n\/[k\(n — k)\]) grows as long as < rrij < K — 1, but decreases again when 
K < mj < 1{K — 1). As a consequence, for any observable there exists a maximal 
m max SUCu that contributions from paths with rrij > mf iax < 2(K — 1) could safely 
be disregarded in the numerical iteration. Of course m'f iax is chosen depending on the 
model parameters and the observable under investigation. 

Rapidly decreasing weights of the paths may not be (over-) compensated by 
increasing weights for < rrij < K — 1, since each contribution is small and the number 
of paths decreases again for larger rrij > K. The behaviour of the impurity weights is 
illustrated as follows. Consider the case when rrij is close to the maximum 2(K — f). 
Both path classes with rrij = and rrij = 2(K — I) contain the same number of elements 
(four), while each path contribution in the second class is weighted by (J5t/2) 2<yK ~ l \ , 
For typical values of K = 4, 5tT = 1/2, and J — T, the weight is ~ 2.5 x I0~ 4 . This 
also holds for all K < rrij < 2(K — I). Since m™ ax is unknown a priori, we include it 
into our code as an additional parameter. Then, we perform a numerical estimate by 
a spot sample of the parameter space. It turns out that for many cases, it is sufficient 
already to choose m™ ax = 2. Fig. 10 shows an example where both (/) and (r z ) converge 
quickly for increasing m™ ax . This drastically reduces the CPU running times, e.g., for 
parameters chosen as in Fig. 10, from more than one month to typically three to five 
days. 
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Figure 11. The expectation value (r z ) (log scale) of the impurity orientation as 
a function of time for four different strengths J of an anti-ferromagnetic electron- 
impurity interaction. The initial preparation of the system at t — — oo is spin-up 
[r z (0) = n = 1] and we set <$> D = A = A imp = U = 0, /3T = 1, and eV = 0.6T. 
The calculated impurity orientation (plot marks) is fitted with good accuracy (errors 
indicated by shaded areas) to exponentially decaying functions (solid lines). 



4. Charge Current and Impurity Dynamics 

The ISPI scheme has originally been developed for the Coulomb-interacting single- 
level quantum dot (SLQD). By reproducing established analytical and experimentally 
confirmed results, the general validity of the approach has been shown in Refs. [48, 49]. 
Here, we focus on novel transport features caused by the magnetic impurity and its 
interaction with dot electrons. We emphasize that novel dynamical and transport 
features are mediated by the transverse or flip-flop interaction H^t, given by Eq. (4). 
Without the possibility of flip-flops the orientation of the impurity spin and its quantum 
state could not change and not participate in the dynamics. The remaining longitudinal 
part of the interaction Hf nt causes a renormalisation of rates and energies which adds to 
the effect stemming from a magnetic field. Necessarily, flip-flop processes are involved 
from the beginning to investigate the non-trivial impurity dynamics by considering the 
time dependence of the impurity orientation (t z ). 

In all presented results below, the impurity is initially polarized and then the 
coupling to the leads is switched on. We observe that the time-decay of the polarisation 
is well described by a single exponential with a constant decay rate. In order to single out 
the relevant physical processes, we compare our numerical ISPI results to a diagrammatic 
perturbation theory in the weak- to intermediate exchange interaction regime. 

We show that for the appropriate parameter regime, the exact numerical results 
are in accordance with the perturbative result and, by this, we obtain a first intuitive 
explanation of the impurity dynamics. A next step is the transfer of its plausibility to the 
ISPI results. However, interaction-induced deviations from the perturbative theory are 
large enough to clearly illustrate the need for a non-perturbative theoretical description, 
provided by the ISPI results. 
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4-1- Real-time decay of the impurity polarisation 

In Fig. 11 we present the time evolution of the impurity polarisation (r 2 ) for different 
values of the exchange interaction J. The remaining model parameters are chosen 
as = A = A imp = U — 0, and (5T = 1 as well as eV = O.QT. As a function 
of time, the impurity polarisation shows a decaying behaviour, well described by an 
exponential relaxation for intermediate to long propagation times. A faster decay of 
the polarisation is observed as the impurity interacts stronger with the electron spins. 
In passing, we note that a rate equation ansatz with constant rate r^ 1 , where tr is 
the relaxation time, results in an exponential decay as well, i.e., (r z )(t) oc e^^^R . 
The parameters are chosen to yield an isotropic (symmetric with respect to [relative] 
spin orientations) model system. In this case the anti-ferromagnetic interaction favours 
anti-parallel orientation of electron- and impurity spin. Over long propagation times, 
the coupling to the unpolarised leads then destroys any polarisation of the impurity. 
It is therefore reasonable to assume, that the rates for up- and down flips are equal. 
In the chosen parameter range, the polarisation of the impurity in contact with the 
leads is well described by a Markovian dynamics, i.e., solely by the time dependent 
probabilities P T (£) of finding the impurity in state |r) at time t. Apparently, this simple 
theoretical prediction agrees well with the numerical results, see Fig. 11. While the 
impurity interaction energy is comparable to the tunnel coupling and considerably affects 
the transport behaviour as we show below (see Fig. 13), the rather high temperature 
and bias voltage nevertheless reduce the relevance of coherent dynamics due to on-dot 
interactions to a secondary role. 

We then turn to the investigation of the inverse relaxation time r^ 1 . In Fig. 12(a), 
we present results for varying J and U — 0, and three different bias voltages. These 
show a nearly quadratic behaviour growing from zero (no relaxation) in the sense that 
for a fit of the results for 0<J<T/2toa polynomial function aJ b the exponent b lies 
between ~ 1.8 and ~ 1.9. An exact quadratic dependence of r^ 1 on J is obtained only 
in cases where the dynamics is strongly dominated by sequential (incoherent) flip-flop 
processes. This is only realized when J <C T. The corresponding rates may be obtained 
based on the real-time diagrammatic technique developed by Schoeller and Schon [24] 
yielding 

r -i^Z!E!v fl, [/l + h+/r + H][/lM + / r hi ( ^ 

Details of the derivation can be found in Ref . [53] . It reveals the physical structure and 
allows for the intuitive interpretation of the processes contributing to sequential flip- 
flops. In the numerator of the integrand, we have the sum of all four possible ways to 
multiply one of the lead's occupations (/ + ) with another or the same lead's probability 
to find an empty state (/ _ ) at some energy. Each of these four combinations is then 
multiplied by the Lorentzian spectral density for the two different spin states each shifted 
by ±J (longitudinal interaction energy). This suggests the following interpretation: A 
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Figure 12. (a) Impurity relaxation rate t^ 1 versus interaction strength J, for three 
different values of V. The parameters correspond to those from Fig. 11. The solid 
lines are fits to polynomial functions aJ b and all resulting values of b are close to 2 
(1.8 < b < 1.9). The polarisation decays faster with increasing J. (b) Comparison 
of the numerically exact [ISPI, crosses] and the sequential relaxation rate [Eq. (54), 
blue solid line] versus bias voltage for J — T/2 and j3T — 1 (here and in what follows, 
parameters not explicitly given are set to zero). Quantitatively, the sequential and 
ISPI relaxation times differ noticeably in the crossover regime. 

sequential flip-flop process consists of three elementary components: the actual flip-flop 
and two tunnelling processes of single electrons with opposite spin (not necessarily in 
that order). Since they evolve coherently, these components form an effective spin-flip 
process \Xi T ) \Xi~ T )i where \ £ {0,cr, d} and the underlying flip-flop nature is 
masked by the tunnelling electrons. For a particular choice of a and a' in Eq. (54), we 
assign certain effective flip processes. 

Fig. 12(b) shows, how the ISPI result (blue crosses) compares to the sequential 
relaxation time (blue solid line). Although the latter is of the correct order of magnitude, 
it is systematically larger than the exact value by > 10%. Since J is not a small 
parameter of the system, we can presume that the deviations are mostly coming from 
coherent higher-order flip-flop processes, which are neglected in Eq. (54). Another source 
of those deviations may be that free Green's functions are used for the derivation of the 
rates in Eq. (54). In their qualitative features, however, both results agree. From their 
finite value at zero bias voltage they grow monotonically. While for small voltages the 
relaxation shows a nearly quadratic functional form (power-law), for larger bias voltages 
eV > 1.25F", it exhibits a linear behaviour. 

4-2. Impact of the impurity interaction on the current 

As opposed to the slow impurity dynamics, measured in terms of the current is 
relaxing fast into the stationary state. This behaviour is caused by the strong coupling 
between the leads and the dot. For the parameters considered here, the upper limit for 
reaching stationarity is about t ST < T^ 1 . Therefore, we consider the stationary current. 
Fig. 13 depicts the current as a function of J with V = 0.6T and (3T = 1 both for 
vanishing Coulomb interaction and for U = T/2. The current decreases with stronger 
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Figure 13. Stationary charge current (I) in units of Iq — eT/h at eV = 0.6L against 
J for two values of U. For increasing interaction, both the (U = 0) Landauer-Biittikcr 
theory and the numerics predict a decrease of the current in the range < J < P. 
The similar characteristics of the LB curve and the ISPI data points suggest that the 
current is mainly affected by the longitudinal part of the electron-impurity interaction. 
This is probably due to the relatively high temperature and, consequently, a short 
coherence time, which strongly limits the influence of coherent dynamics. While the 
LB theory and ISPI agree for J = L/8, growing differences for increasing J show the 
effect of flip-flop processes. The current for small finite Coulomb interaction (no error 
bars given), though consistently smaller than the LB values and also decreasing with 
J, drops slower than for vanishing U. 



impurity interaction. To distinguish the influence of the longitudinal (single-particle) 
and transversal (spin-scattering) part of the interaction, we compare the ISPI results 
with the Landauer-Biittiker (LB) current (/)lb (see [56]), which can be written here as 



For J ^ this is an approximate expression, as it only includes the effect of the 
longitudinal impurity interaction, which acts as an effective magnetic field in the sense of 
a mean field. Similar to the sequential relaxation rates of Eq. (54), the current formula 
has a simple physical interpretation. The joint density of dot-electron states is given by 
a Breit-Wigner function, whose width equals the tunnel coupling strength and whose 
resonance lies at the single-electron energy u a ± J. Hence, the (non-interacting) current 
is given by the integral over the energy-dependent difference of the left and right lead's 
occupation multiplied by the density of available dot states at that same energy. The 
difference in occupation of the lead electronic states is largest around the Fermi level, 
where it has the biggest overlap with the density of states for J = 0. With increasing J, 
the density resonances "move away" from the Fermi level, where f^(u)) — f^(u) decreases 
and the current drops. This effect is explained in terms of the single-particle energy shift 
due to the longitudinal component of the impurity interaction only. 

Fig. 13 shows that the flip-flop term H^ t has a considerably smaller influence on the 
charge current at this rather large temperature (incoherent regime) than the longitudinal 
part of the interaction. Despite the qualitatively similar behaviour of the LB current and 
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Figure 14. Relaxation rate versus bias voltage for three different temperatures and 
U = r/2, J = r. Shown are the ISPI data (symbols) and the perturbative (PT) results 
(solid lines). When compared to the results of Fig. 12(b), where [7 = and smaller J, 
we see larger relative deviations between the ISPI results and the sequential flip-flop 
approximation. The qualitative differences are also more pronounced. 

the exact data, the flip-flop scattering causes an additional significant current drop that 
grows for growing J. A finite Coulomb interaction of U = r/2 increases the resistivity 
of the dot and the ISPI points are consistently lower than the LB values. The decreasing 
effect of the longitudinal part of the impurity, however, is partially compensated by a 
broadening of the joint density of states due to Coulomb fluctuations. 

4-3. Finite Impurity Interaction and Coulomb Repulsion 

In the deep quantum regime, where no small parameter exists, ISPI is certainly 
applicable and able to describe physical properties not predictable by perturbative 
methods. In this section we study, how the relaxation rate and the current behave 
as functions of bias voltage, Coulomb interaction and temperature, respectively. 

Figure 14 presents results in a voltage range < eV < 3T, with J — T, 
U = r/2 and for temperatures 0F = 1, 2 and 5. The ISPI data of r^ 1 are 
indicated by the symbols, while the solid lines mark the corresponding perturbative 
rates. The latter exhibit the same features (power-law growth, followed by a [quasi- 
linear behaviour, which finally saturates) as in Fig. 12(b), which are more pronounced 
for lower temperatures. As expected, the ISPI data points deviate considerably from 
this lowest-order approximation. Quantum coherent effects are increasingly relevant 
since, all energies are of the same order. Both, the degree of quantitative differences 
and the deviations in the qualitative behaviour increase with lower temperatures. 

In Fig. 15, for each of two different temperatures four different current curves are 
shown — one for each possibility to either have (i) only mean field dynamics (LB), (ii) the 
full Coulomb interaction without flip-flop processes ("no flips"), (iii) flip-flop dynamics 
without Coulomb fluctuations ("mean-field U"), or (iv) the fully interacting dot ("full 
int."). For J = T and V = 2T, the Coulomb energy is varied between < U < T. The 
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Figure 15. Comparison of (i) the LB current (dashed lines), (ii) the Coulomb 
interacting current without flip-flop scattering ("no flips"), (hi) the current without 
Coulomb scattering but full impurity interaction ("mean-field U", see text), and (iv) 
the fully interacting current ("full int.") in their dependence on the Coulomb interaction 
U for (a) j3T — 1 and (b) (3T = 5. The other (non-zero) parameters are J = T and 
V = 2r. The red-shaded areas indicate the error margin for the case "mean-field £/". 



situation "mean- field U" is implemented by setting $£> = U/2 and the HS-parameter 
A = to illustrate the effect of the "classical" or mean-field part of the Coulomb 
interaction. This is tantamount to setting the third term in Eq. (11) to zero thereby 
neglecting the Coulomb interaction induced fluctuations, which results in a shift of each 
single electron energy by U/2. Since this shift tends to move the transport channels 
away from the Fermi level, i.e., the region with the highest density of state in the leads, 
the current drops. By its nature, this decrease is equivalent to the one observed for the 
LB current. 

Only for the "single-interaction" currents ("no flips"), we show the error bars. The 
reason why no margin of confidence is given for the fully interacting case, regards the 
comparability of the error data. Calculating the "full int." current is a very time 
consuming task and thus, the extrapolation involves considerably fewer data points. 
Nevertheless, this does not render these values unreliable. We still see a compelling 
linear behaviour of the l/r c extrapolation with errors of the order of 1% based on the 
sample standard deviation. Notice that with about 10%, the relative error of the current 
values is rather small, the small variations of the "full int." data are solely due to the 
extrapolation errors and have no physical meaning. 

For both temperatures, both the LB current and the current without Coulomb 
scattering show only a weak dependence on U due to the single-particle energy shift. 
The current with full Coulomb interaction but fixed impurity shows a local maximum 
for U ~ T/2, which is more pronounced for /3T = 5. In this case, the fixed impurity acts 
as an effective static magnetic field. The ISPI values for the fully interacting dot vary 
strongly over the considered U interval, but are scattered around the "no flips" and "no 
U" curves. 

As long as the Coulomb interaction is small, all current values in both figures lie 
close. The difference between the respective current values is given by the inclusion or 
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exclusion of flip-flop processes only. Hence the rather good agreement of the U — values 
suggests, that even at this temperature flip-flop processes alone affect the current only 
weakly. Nevertheless, for lower temperatures, the flip-flop processes start to influence 
the current more strongly, which results in an increased resistivity. The case of the "no 
flip" current (fixed impurity) is equivalent to a Coulomb-interacting SLQD in a magnetic 
field. Both curves in (a) and (b) show a very similar dependence on U, featuring a local 
maximum at around U = T/2. For the lower temperature, the relative height of the 
broad current peak is twice as big as for j3T = 1. This effect is also caused by the 
broadening of the dot's joint density of states due to the Coulomb fluctuations. 

5. Conclusions 

We have studied the real-time nonequilibrium dynamics of a single-orbital magnetic 
quantum dot including Coulomb interactions. In order to obtain stationary 
nonequilibrium states at asymptotic times, the ISPI scheme is employed and extended to 
cases, when an additional magnetic degree of freedom is present. Besides the appearance 
of a Hubbard-Stratonovich field, which decouples the Anderson repulsion term in the 
Hamiltonian, we have to include the impurity interaction on the same level. This 
nontrivial task requires an additional summation over paths of the impurity spin degree 
of freedom. The resulting action in the path integral formalism involves the Green's 
function of the quantum dot as well as its inverse. Inversion of the Green's matrix 
enlarges the numerical effort tremendously. From the technical point of view, appearing 
matrices, dealing with the impurity dynamics may violate the necessary block diagonal 
structure. However, a unitary transformation helps to build up proper ingredients for the 
algorithm. Then, also impurity induced correlations become tractable and do not violate 
the exponential decay of quantum many-body correlations. We have presented how 
an efficient truncation scheme provides accurate results for the coupled spin dynamics. 
Results are given for a quantum spin-1/2 impurity on the dot, whereas the generalization 
to an impurity with a larger spin is possible. This would be necessary when, e.g. a Mn 
system is under investigation. 

For the same kind of exchange coupling, the implementation of a Mn impurity 
essentially increases the dimension of the impurity path sum. Instead of summing over 
all step-like paths of a spin-1/2, it involves the sub-class of step-like paths (steps between 
orientations differing by one) in the space of the six possible orientations of spin 5/2. Of 
course, the numerical efforts also increase but still are within reach of the ISPI method. 
Further work will be dedicated to this goal. 

The exponential drop of time correlations due to the leads' coupling allows for an 
efficient truncation of the appearing sums — the main building block of the ISPI scheme. 
Its application yields high quality numerical data for the impurity relaxation time and 
the tunnelling current as a function of the bias voltage in the presence of (Coulomb) 
and electron-impurity interactions. We have performed necessary checks to compare 
our findings to established results. In the regime of small impurity interaction, where 
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sequential flip-flops dominate the impurity dynamics, we have found good agreement 
with a classical rate equation. This is a useful tool to gain insight into the dominating 
processes in the incoherent regime. Relaxation is described reasonably well by a rate 
equation when lead-induced coherences are absent. 

In the deep quantum regime, however, we find that the ISPI method is the only 
tool to obtain both the correct order of magnitude and the qualitative features of the 
relaxation rate as it depends on the system parameters U and J. The same holds for 
the influence of J and U on the current in the deep quantum regime. Most importantly, 
the ISPI scheme proves to be useful to cover the full cross-over regime where no small 
parameter can be identified and thus any perturbative approach becomes invalid. 

Furthermore, Kondo physics in such a single spin system under nonequilibrium 
conditions is of course an interesting subject to study. It emerges when the Coulomb 
interaction is large compared to the tunnel coupling and the temperature is sufficiently 
low (also in comparison to the tunnel coupling). In the present work, the two interaction 
terms in the Hamiltonian and the strong tunnel coupling presently limit the application 
of ISPI to intermediate temperatures. Therefore, Kondo features have not yet been 
obtained. In that regard, the further development of the method is still demanding, 
see also the discussion in Ref. [48]. Nevertheless, due to the suppression of long- 
time correlations at finite voltages, the regime of nonlinear transport where the Kondo 
correlations compete with the finite bias is accessible and will be treated in future work. 

We have provided a first glimpse on the interesting new physics that comes into 
reach with the ISPI scheme. A generalisation of the model to several localised magnetic 
impurities, with electrons mediating a finite magnetisation between them should be 
possible. The real-time dynamics and all-electrical control of such devices could be 
simulated. The presented scheme is also applicable to provide the x- and y-components 
of the impurity spin, thus yielding the complete spin dynamics and the real time 
dephasing on the Bloch sphere. 
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